Seismic Inverse Scattering via Discrete Helmholtz Operator Factorization and Optimization

We present a joint seismic inverse scattering and finite-frequency (reflection) tomography program, formulated as a coupled set of optimization problems, in terms of inhomogeneous Helmholtz equations. We use a higher order finite difference scheme for these Helmholtz equations to guarantee sufficient accuracy. We adapt a structured approximate direct solver for the relevant systems of algebraic equations, which addresses storage requirements through compression, to yield a complexity for computing the gradients or images in the optimization problems that consists of two parts, viz., the cost for all the matrix factorizations which is O(rN log N) times the number of frequencies, and the cost for all solutions by substitution which is O(N log(r log N)) times the number of frequencies times the number of sources (events), where N = n if n is the number of grid samples in any direction, and r is a parameter depending on the preset accuracy and the problem at hand. With this complexity, the multi-frequency approach to inverse scattering and finite-frequency tomography becomes computationally feasible with large data sets, in dimensions d = 2 and 3.